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Abstract 

We study the finite density phase transition in the lattice QCD at real chemical potential. We 
adopt canonical approach and the canonical partition function is constructed for Nj = 2 QCD. 
After derivation of the canonical partition function we calculate observables like the pressure, the 
quark number density, its second cumulant and the chiral condensate as a function of the real 
chemical potential. We covered a wide range of temperatnre region starting from the confining 
low to the deconfining high temperature. We observe signals for the deconfinement and the chiral 
restoration phase transition at real chemical potential below Tc starting from the confining phase. 
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I. INTRODUCTION 


The grand canonical ensemble is a difficult subject to treat in lattice QCD because of 
the sign problem. The canonical partition function is related to the grand canonical one 
through the fugacity expansion and is known to be free from the complex action problem. 
The history of the canonical approach in hnite density lattice QCD may have started with 
a development of a reduction formula for the Dirac determinant [l|, which gives naturally a 


fugacity expansion of the Dirac determinant. A several studies with the staggered 


are done along this line 


^ and have been taken over by the Wilson fermion 


10 


ermion 
. Since 


it is understood that the inverse transformation of the fug acity expansion is the Fourier 


transformation 


3j it also becomes a popular method 
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id- 


One of the problem of the canonical approach is that the Dirac determinant needs to be 
evaluated accurately, which is the heaviest part of the computation. One of the solution 


is to adopt the Taylor expansion 1^. In this paper we perform the fugacity expansion 
by a method of the hopping parameter expansion in temporal direction: winding number 
expansion j^, 13|. By using the hopping parameter expansion the Dirac determinant can be 
evaluated with low cost and we can visit a wide parameter space in f3. The canonical partition 
function is evaluated for Nf = 2 QCD both in the deconhnement and the confinement 
temperature regions. After derivation of the canonical partition function we study the 
chemical potential dependence of observables like the pressure, the quark number density 
and the chiral condensate. We observe signals for the confinement-deconfinement and chiral 


restoration phase transition at real chemical po 


preliminary result has been reported in Ref. 


li 


ential below the critical temperature T^. A 


This paper is organized as follows. In Sec. 2 we briefly review the canonical approach. 
Our formulation of the fugacity expansion in terms of the hopping parameter expansion is 
discussed in Sec. 3. After mentioning our numerical setup in Sec. 4 our main results are 
given in Sec. 5 for the canonical and grand canonical partition functions. Sec. 6 is devoted 
for the chiral condensate in grand canonical ensemble. A conclusion is given in Sec. 7. 
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II. CANONICAL PARTITION FUNCTION 


It is a well known fact that the grand canonical ensemble and the canonical one are 
equivalent each other. This is shown by a simple equation to relate the grand canonical 
partition function T, V) and the canonical Zc{n,T, V) 

OO 

ZG{^^,T,V)= Zc{n,T,V)C, i = ( 1 ) 


This is a so called fugacity expansion formula. The inverse of this expansion is given by 
using the Cauchy’s integral theorem 


Zo (n,T,V) = j, T, V), 


( 2 ) 


where we assume that the partition function Zg{^,T, V) has singularities only at .^ = 0, cxo 
corresponding to trivial ones /i/T = ±oo and adopt an appropriate contour around the 
origin. Here we notice that a phase transition point is not a singularity of the partition 
function but it is rather a zero of ZqUiT, V) (Lee-Yang zeros 1^, ll9|) and does not affect 
the Cauchy’s integral. 

Now it is free to change the contour to a unit circle ^ = e*® and the contour integral turns 
out to be a Fourier transformation [Si 


f2TT fja 

Zc{n,T,V)= / —e-'"®ZG(e*^T,U). 

Jo 2 tt 


( 3 ) 


A standard Monte Carlo simulation is possible for the lattice QCD since the chemical po¬ 
tential is set to pure imaginary /r/T = i6 and the Boltzmann weight is real positive for two 
flavors. 

However a difficulty of the sign problem is preserved unfortunately since there remains a 
frequent cancellation in plus and minus sign in a numerical Fourier transformation especially 
for large particle numbers. In order to get the canonical partition function accurately for 
large quark number n we need a very hne sampling of Zoie^^,T,V) in 6. This requires a 
heavy computational cost because we need to evaluate the Dirac operator determinant for 
the QCD grand partition function. In this paper we shall solve this problem by a direct 
expansion of the Dirac determinant in terms of the fugacity. 
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III. WINDING NUMBER EXPANSION 


We consider the lattice Nf = 2 QCD grand partition function given in the path integral 
form 


Zg{^,T,V) = J 


(4) 


where we adopt the improved Wilson Dirac operator 
D\y = 1 — nQ, 

Q = E {Qt + Qr) + + e-i^Q^ + T, 


i=l 


— (1 7Ai) U^{n)5m,n+[ii ~ 4" 7/^) Ujj^(jn)6rn,n-f},, 

i'^)nm ~ '"SWTT F^^(rL)5 mn 


(5) 

( 6 ) 

(7) 

( 8 ) 


ll,U 


Both the chemical potential and the hopping parameter k appear in front of the 
temporal hopping term simultaneously. The fugacity expansion of the determinant shall be 
given by using the hopping parameter expansion. 

Instead of the determinant we consider the hopping parameter expansion of 


oo 


TrLogDvi/ = TrLog (/ — kQ) = — ^ —Tr(Q 


n=l 


n 


(9) 


A non-trivial chemical potential dependence appears when the quark hoppings make a loop 
in the temporal direction. If one of the term has a n times winding loop in positive temporal 
direction, the chemical potential dependence is where Nt is a temporal lattice length 

and this is nothing but Counting the winding number in temporal direction for 

each term in equation ([9]), we get the winding number expansion of the TrLog 

OO 

TrLogT>vc(/i) = WnC- (10) 

n=—oo 

Regrouping the summation, we have a fugacity expansion of the determinant [l2, ^ 

/ OO \ OO 

DetDH.(e;t/) = exp ^ WnC] = E ^niU)C. (11) 


In this approach, the fugacity dependence of the determinant or the partition function is 
given analytically. A numerical Fourier transformation ([3]) can be executed safely at any 
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/S 

K, 

m-„/mp 

0.9 

0.137 

0.8978(55) 

1.1 

0.133 

0.9038(56) 

1.3 

0.138 

0.809(12) 

1.5 

0.136 

0.756(13) 

1.7 

0.129 

0.770(13) 

1.9 

0.125 

0.714(15) 

2.1 

0.122 

0.836(47) 


TABLE I: (3 and k for the numerical simulation. rriT^/mp is also given. 

high precision we want. In this paper, extraction of the coefficient Zn{U) from the second 
term in eqnation ffTT]) is done by using the multi-precision numerical Fourier transformation. 

The winding number expansion is done for gauge conhgurations generated at /xq = 0 or 
purely imaginary value. This procedure corresponds to the standard reweighting method 

z,K.r. V) = I DU ( 12 ) 

IV. NUMERICAL SETUP 


We adopt the Iwasaki gauge action and the improved Wilson fermion action with the 
clover term. The number of flavors is set to Nj = 2 with degenerate masses. The APE stout 


smeared gauge link 


20| is used for those in the fermion action including the clover term. 


The number of smearing is four and the parameter is set to p = 0.1. The clover coefficient 
is hxed to csw = 1-1 (2l|. We adopt 8^ x 4 lattice. A wide range of {3 is covered from high 
temperature jS = 2.1 to low temperature 0.9. It seems to be that both the conhning and 
deconhning region are well covered, which is inferred from a behavior of real part of the 
Polyakov loop in the left panel of Fig. [1] The temperature corresponding to /3 = 1.7 seems 
to be very near to the critical as is shown the right panel of Fig. [1] using a fluctuation of 
phase of the Polyakov loop. The hopping parameter is selected in order that the hopping 
parameter expansion works well, for which m.,^/mp turns out to be 0.7 - 0.9 as is given in 
table [H Maximal order of the hopping parameter expansion is set to 480 so that max winding 
number in temporal direction is 120. Number of independent conhguration is 100 - 600. 
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FIG. 1: Real part of the Polyakov loop as a function of /3 (left). Both the confining and deconfining 
regions are well covered. Right panel is a fluctuation of phase of the Polyakov loop. /3 = 1.7 seems 
to be very near to the critical coupling. 
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FIG. 2; The canonical partition function log |Zc(nB)/.^c(0)|/(RT^) as a function of the baryon 
number ub- From the top j3 = 2.1 (orange), 1.9 (red), 1.7 (magenta), 1.5 (green), 1.3 (dark green), 
1.1 (blue), 0.9 (black). 


V. CANONICAL AND GRAND CANONICAL PARTITION FUNCTION 


The first numerical result we get is the canonical partition function Zc{n^ T, V). We plot 
log \Zc{nB-, T, l/)/Zc(0, T, V)\/{VT^) as a function of the baryon number ub in Fig. [21 The 
partition function decays very rapidly with ub and its behavior changes drastically between 
(3 = 1.9 (red) and 1.5 (green), which may correspond to a phase transition. 

The plot range is hxed by using the d’Alembert’s convergence condition 


lim 

riB^OQ 


Zcij^B + 1 ) 

Zc{nB) 


(13) 
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FIG. 3: Log plot of \Zc{nB + 1)/Zc{nB)\ as a function of the baryon number. The color and P 
correspondence is the same as in Fig. 2 for left panel. Right panel is that at /3 = 1.1 and 2.1 for 
example. 


which gives the convergence radius for the fugacity ^ The data at (3 = 1.1 and 2.1 are 
plotted in right panel of Fig. [3] for example. We cut our data at Umax where a monotonic 
decrease stops indicated by vertical orange and blue lines in the figure. The horizontal lines 
show maximal values of the fugacity expected to be within the convergence radius. By 
taking log, the line gives our applicable limit for the baryon chemical potential /is/T. For 
example we can discuss physics safely at —10 < ps/T < 10 for f3 = 1.1 and —4 < Hb/T < 4 
for (3 = 2.1 with our method. We notice that Umax is improved only by one even if we 
increase the maximal winding number twice as 120 —?■ 240. This may be due to a fact that 
the convergence radius of the Taylor expansion ([3]) is unity. 

The second physical quantity is the grand partition function. By taking summation for 
—'^ma.x ^^B < Umax in cquatiou ([T]), wc get the grand partition function for the real chemical 
potential. Here we notice that the canonical partition function Zc{n, T, V) is a real positive 
quantity if a Hermitian transfer matrix exists for the lattice QCD. We observed that the 
absolute value of the partition function is very stable against the statistical fluctuation and 
we consider it is suitable for extracting physical information. Instead of equation ([T]) we 
adopt 

^max 

Zg(m,T,V)= ^ \Zc{n,T,V)\C (14) 

^max 

^ For negative baryon number we adopt \Zc{nB — l)/ZcinB)\. 
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FIG. 4: The grand partition function log T, I/)/Zg'( 0, T, I/)|/(FT^) as a function of the 

quark chemical potential ^/T. The color and /3 correspondence is the same as in Fig. 2. The 
pressure at /3 = 0.9 and 1.1 are plotted in the upper left panel, where two data are completely 
overlapped. The upper right panel is that at /3 = 1.3 and 1.5. Lower panel is that for /3 = 1.7, 1.9, 

2 . 1 . 


for a definition of the grand partition fnnction from which we get the pressnre in the grand 
canonical ensemble 


T(P(^/r)-P(0)) = pGiog 


Zcip,T,V) 


(15) 


Zo(0,T,V) 

We plot the pressnre normalized at /i = 0 in Fig. H] as a fnnction of the qnark chemical 
potential. We observe that the pressnre is almost consistent with zero for small chemical 
potential at low temperatnre below Tc as is shown in the npper panels of Fig. IH 

The third physical application is the k-th momentnm of the qnark nnmber operator 

1 


(N^) = 


Zg{p,T,V)^Jz 


Y. n^Zcin,T,V)e 


(16) 


Notice this is not a phase quenching of the Dirac determinant. 
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FIG. 5: Quark number density (N)/{VT^) as a function of the quark chemical potential. The 
color and /3 correspondence is the same as in Fig. 2. The upper panels plot results for /3 = 0.9 - 
1.5 below Tc- The lower panel is results for f3 = 1.7 - 2.1 above Tc. 


In Fig. [5l the quark number density {N) / (VT^) is plotted as a function of the real quark 
chemical potential p/T. We observe a Silver blaze like phenomenon at low temperature 
below Tc as is shown in the upper panels. The quark number density is consistent with 
zero within a statistical error until p/T ~ 1 is reached. Above p/T ~ 1 the quark number 
density takes off and becomes non-zero, which indicates deconhnement phase transition. 

The second cumulant of the quark number density may be a good indicator of this phase 
transition. The results at low temperature side below Tc are plotted in the upper panels of 
Fig. [6l Those above Tc are given in the lower panel. We observe peaks around p/T ~ 1 
below Tc for every (3 in the upper panels. We also observe secondary peaks for (3 = 1.3 
(dark-green) and 1.5 (green) at larger p/T > 2 in the upper right panel. Only a single peak 
is observed at high temperature in the lower panel. These peaks have a distinctive behavior 
against the quark number cut-off nmax- The hrst peaks below Tc are stable even if we vary 
^max by 20%. On the other hand, the second peaks at (3 = 1.3, 1.5 and the hrst peaks above 
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FIG. 6: The second cumulant of the quark number (iV^)c/ iVT^) as a function of the quark chemical 
potential. The upper panels plot results for /? = 0.9 - 1.5 below Tc- The lower panel is results for 
/3 = 1.7 - 2.1 above Tc. 


Tc change their height and position with rimax- It would be conjectured that the hrst peaks 
in the low temperature side are the physical indication of the conhnement-deconhnement 
phase transition. The secondary peaks below Tc and the peaks in the high temperature side 
would be artifacts due to the cut-off in quark number in the fugacity expansion. 

The thermodynamical quantity 

/i d\nZc{n,T,V) I d\nZc{pV,T,V) 


T 


On 


1 / 


dp 


(17) 


as a function of the 
canonical approach 


B 


jaryon number density is one of the post popular quantity in the 


13 


16| . We notice the quantity plotted in Fig. [3] can be interpreted 


as an expression of the derivative in terms of the differential 

I = -1 (In Zc(nB + 1, r, C) - In Zc(nB, T,V)). 


(18) 


The results are plotted in Fig. [3 Although the error is rather large there seems to appear 
a typical S-shape like behavior at small baryon density Pb/T < 0.5 for the low temperature 
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FIG. 7: The quark chemical potential ^/T given in terms of the thermodynamical relation as a 
function of the baryon number density. The upper panels plot results for (3 = 0.9 - 1.5 below T^. 
The lower panel is results for ^0 = 1.7 - 2.1 above Tc. 


side given in the upper panels, which may indicate an existence of the hrst order phase 
transition. 


VI. HADRONIC OBSERVABLES 

In our procedure, the fugacity expansion is based on the hopping parameter expansion. 
It may be possible to expand any hadronic operators in terms of the fugacity. We consider 
a fugacity expansion of a numerator of some operator vacuum expectation value (VEV) 

/ OO 

Oc(n,T,V)e. (20) 

n=—oo 

where (O)quark ^ Wick contraction of a hadronic operator O in quark helds. For example, 
we consider the chiral condensate. It is easy to expand its Wick contraction in terms of the 
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hopping parameter 

(V’V') = “tr (-Fr~) = 1 I 

\ / quark V J \ 1 “ / 


^«:-trg-= ^ o.(f/)r. 

m=0 n=—oo 


( 21 ) 


Counting the winding number in temporal direction, we get the last equality. Multiplying 
the determinant contribution (ITT]) , we apply the same regrouping procedure in Sec. IIIII and 
get the fugacity expansion fl20l) . 

Once we have two coefficients Zc{n,T,V) and Oc{n,T,V), the VEVs of the operator 
with the canonical and the grand canonical ensemble are available. The canonical ensemble 
VEV is given by taking the ratio of two coefficient 


(0)c (n, r, V) 


Oc(n,T,V) 

Zc(n,T,V)' 


( 22 ) 


The result for the chiral condensate — / d^x / (VT^) is given in the left panel of Fig. [8] 

as a function of the baryon number. 

A VEV in the grand canonical ensemble is given by taking fugacity summation with real 
chemical potential 


(0)^(e,T,V) = 




(23) 


The chiral condensate — f d^x / (VT^) is given in the right panel of Fig. |8]as a function 

of the real chemical potential p/T. The condensate in the hgure is a bare quantity without 
renormalization. Since we adopted the Wilson fermion, we have an additive correction for 
which is not subtracted in this paper. From the hgure, the chiral restoration phase 
transition at hnite chemical potential seems to be seen. A relatively large value at around 
p/T = 0 starts to decrease rapidly at chemical potential p/T ~ 1 for low temperature. The 
would-be transition parameter Pc/T becomes larger for lower temperature and smaller for 
high temperature although we need to notice the quark mass is not the same for each fd. The 
position of ^c/T seems to be consistent with that of the hrst peak for the second cumulant 
in the upper panels of Fig. |6l 

The quark number density is also given in terms of the hadronic observable 
/ dd" /(VT^), which gives consistent result with that in Sec. 5 up to renormalization 

factor as is plotted in FigO 
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FIG. 8: The chiral condensate — / d^x{ipip) /{VT^) in the canonical ensemble (left panel) as a 
function of the baryon number. The right panel is that in the grand canonical ensemble as a 


function of the quark chemical potential /x/T. The color and f3 correspondence is the same as in 
Fig. 2. 



FIG. 9: The quark number density / d^x /iVT^) in the grand canonical ensemble as a 

function of the quark chemical potential ^/T. The color and /? correspondence is the same as in 
Fig. 2. 


VII. CONCLUSION 

In this paper, we performed the fugacity expansion of the grand partition fnnction by 
using the hopping parameter expansion. This procedure seems to be valid for baryon num¬ 
bers around ub ~ 30 for 8^ x 4 lattice and m-,^/mp > 0.7. The method is also applied to the 
numerator of VEVs of hadronic operators. Taking summation over quark numbers we get 
a VEV at the real chemical potential. As an example we evaluate the pressure, the quark 
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number density and its second cumulant and the chiral condensate. These observables show 
a phase transition like behavior at high chemical potential for low temperature region, which 
may correspond to the conhnement-deconhnement or the chiral restoration phase transition. 

One of the biggest anxiety in this paper may be an overlap problem due to the reweighting 
m- We performed a test of this problem by adopting po/T = 0.5 for conhgurations 
generation. We hnd all the results are consistent within the statistical error as we shall 
report in the forthcoming paper. 
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